## script to reclassify altitude output rasters to desired values and to create a flat lanscape
# plot rasters afterwards

rm(list = ls())
library(sp)
library(raster)
library(rgdal)
## rgdal: version: 1.5-10, (SVN revision 1006)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.4.0, released 2018/12/14
## Path to GDAL shared files: /usr/share/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
## Path to PROJ shared files: (autodetected)
## Linking to sp version:1.4-2
##### Create "flat" landscape #####
source("create_Flat.R")

create_Flat(inraster="../spatial/Elevation/NevTol_Alt.tif",
            outname="../spatial/Elevation/NevTol_Alt_flat.asc")
## [1] "original raster"
## class      : RasterLayer 
## dimensions : 3399, 3582, 12175218  (nrow, ncol, ncell)
## resolution : 9.316244e-05, 9.316384e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /LUSTRE/Genetica/ngalvez/Multihierarchical_NevadoToluca/spatial/Elevation/NevTol_Alt.tif 
## names      : NevTol_Alt 
## 
## [1] "output raster"
## class      : RasterLayer 
## dimensions : 3399, 3582, 12175218  (nrow, ncol, ncell)
## resolution : 9.316244e-05, 9.316244e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : /LUSTRE/Genetica/ngalvez/Multihierarchical_NevadoToluca/spatial/Elevation/NevTol_Alt_flat.asc 
## names      : NevTol_Alt_flat
##### Reclassify models  #####
# so that non suitable habitat is set to NotS and suitable habitat to 1
source("reclass4circuit.R")
datafolder="../spatial/SDM/out"


##### Reclassify elevation models  #####
# so that not suitable elevation is set to NotS and suitable to 1

# path to elevation raw data
datafolder="../spatial/Elevation"

# define in data and function
inraster=paste0(datafolder, "/NevTol_Alt.tif")
source("reclassElev4circuit.R")

## Do reclassification for several elevations

#for (i in c(2000, 2700, 2800, 2900, 3000, 3100, 3200, 3300, 3400, 3500, 3600, 3700)) {
#    print(paste("for min elevation:", i))
#    outraster=paste0(datafolder, "/NevTol_Alt_", i, "_reclass.asc")
#    # reclassify
#    reclassElev4circuit(inraster, minElev=i, outraster, NotS=0.2)
#    }


##### Raster Elevation Hipotesis "A"######### 
myraster<-raster("../spatial/Elevation/NevTol_Alt.tif")
myraster
## class      : RasterLayer 
## dimensions : 3399, 3582, 12175218  (nrow, ncol, ncell)
## resolution : 9.316244e-05, 9.316384e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /LUSTRE/Genetica/ngalvez/Multihierarchical_NevadoToluca/spatial/Elevation/NevTol_Alt.tif 
## names      : NevTol_Alt
plot(myraster)

# read reclasification matrix (make sure no header)
# reclassify the values into three groups 
rcl_Alt_A_reclass <- read.table("../spatial/Elevation/RasterAltitude_A.txt", sep = ",", header=F)
dim(rcl_Alt_A_reclass)
## [1] 9 3
class(rcl_Alt_A_reclass)
## [1] "data.frame"
rcl_Alt_A_reclass<- as.matrix(rcl_Alt_A_reclass) 
class(rcl_Alt_A_reclass)
## [1] "matrix"
rcl_Alt_A_reclass
##         V1   V2  V3
##  [1,] 4000 4457 0.0
##  [2,] 3386 4000 0.1
##  [3,] 3189 3386 0.3
##  [4,] 3164 3189 1.0
##  [5,] 3139 3164 1.0
##  [6,] 3010 3139 0.3
##  [7,] 2600 3010 0.2
##  [8,] 2000 2600 0.1
##  [9,]    0 2000 0.0
# reclasify input raster 
xa<-reclassify(myraster,  rcl=rcl_Alt_A_reclass, include.lowest=FALSE, right=FALSE)
xa
## class      : RasterLayer 
## dimensions : 3399, 3582, 12175218  (nrow, ncol, ncell)
## resolution : 9.316244e-05, 9.316384e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : NevTol_Alt 
## values     : 0, 1  (min, max)
plot(xa)

plot(xa, col=c("grey", "black"), legend=FALSE, xaxt='n', yaxt='n')

# save new raster
writeRaster(xa, filename="../spatial/Elevation/rcl_Alt_A.tif", format="GTiff", overwrite=TRUE)
writeRaster(xa, filename="../spatial/Elevation/rcl_Alt_A.asc", format="ascii", overwrite=TRUE)


###### Raster Altitude Hipotesis "B"############ 
# reclassify the values into three groups 

rcl_Alt_B_reclass <- read.table("../spatial/Elevation/RasterAltitude_B.txt",sep = ",", header=F)
dim(rcl_Alt_B_reclass)
## [1] 8 3
class(rcl_Alt_B_reclass)
## [1] "data.frame"
rcl_Alt_B_reclass<- as.matrix(rcl_Alt_B_reclass) 
class(rcl_Alt_B_reclass)
## [1] "matrix"
rcl_Alt_B_reclass
##        V1   V2  V3
## [1,] 4000 4457 0.0
## [2,] 3386 4000 0.2
## [3,] 3189 3386 0.3
## [4,] 3164 3189 1.0
## [5,] 3139 3164 1.0
## [6,] 3010 3139 0.3
## [7,] 2000 3010 0.2
## [8,]    0 2000 0.0
# reclasify input raster 
xb<-reclassify(myraster,  rcl=rcl_Alt_B_reclass, include.lowest=FALSE, right=FALSE)
xb
## class      : RasterLayer 
## dimensions : 3399, 3582, 12175218  (nrow, ncol, ncell)
## resolution : 9.316244e-05, 9.316384e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : NevTol_Alt 
## values     : 0, 1  (min, max)
plot(xb)

plot(xb, col=c("grey", "black"), legend=FALSE, xaxt='n', yaxt='n')

# save new raster
writeRaster(xb, filename="../spatial/Elevation/rcl_Alt_B.tif", format="GTiff", overwrite=TRUE)
writeRaster(xb, filename="../spatial/Elevation/rcl_Alt_B.asc", format="ascii", overwrite=TRUE)


##################SLOPE#############################


#Raster Slope Nevado de Toluca
myrasterS<-raster("../spatial/Slope/NevTol_Pen.tif")
myrasterS
## class      : RasterLayer 
## dimensions : 3399, 3582, 12175218  (nrow, ncol, ncell)
## resolution : 9.316244e-05, 9.316384e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /LUSTRE/Genetica/ngalvez/Multihierarchical_NevadoToluca/spatial/Slope/NevTol_Pen.tif 
## names      : NevTol_Pen
plot(myrasterS)

#######Raster Slope hipotesis "A"#####################
rcl_S_A <- read.table("../spatial/Slope/RasterValue_SlopeA.txt",sep = ",", header=F)
dim(rcl_S_A)
## [1] 3 3
class(rcl_S_A)
## [1] "data.frame"
rcl_S_A<- as.matrix(rcl_S_A) 
class(rcl_S_A)
## [1] "matrix"
rcl_S_A
##         V1      V2  V3
## [1,]  0.00  6.5700 1.0
## [2,]  6.57 13.1400 0.5
## [3,] 13.14 19.6999 0.2
# reclasify input raster Slope NT
x_SA<-reclassify(myrasterS,  rcl=rcl_S_A, include.lowest=FALSE, right=FALSE)
x_SA
## class      : RasterLayer 
## dimensions : 3399, 3582, 12175218  (nrow, ncol, ncell)
## resolution : 9.316244e-05, 9.316384e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : NevTol_Pen 
## values     : 0.2, 1  (min, max)
plot(x_SA)

# save new raster
writeRaster(x_SA, filename="../spatial/Slope/rcl_S_A.tif", format="GTiff", overwrite=TRUE)
writeRaster(x_SA, filename="../spatial/Slope/rcl_S_A.asc", format="ascii", overwrite=TRUE)


#############Raster Slope Hipotesis "B"###################
rcl_S_B <- read.table("../spatial/Slope/RasterValue_SlopeB.txt",sep = ",", header=F)
dim(rcl_S_B)
## [1] 3 3
class(rcl_S_B)
## [1] "data.frame"
rcl_S_A<- as.matrix(rcl_S_B) 
class(rcl_S_B)
## [1] "data.frame"
rcl_S_B
##      V1      V2  V3
## 1  0.00  6.5700 1.0
## 2  6.57 13.1400 0.7
## 3 13.14 19.6999 0.2
# reclasify input raster Slope NT
x_SB<-reclassify(myrasterS,  rcl=rcl_S_B, include.lowest=FALSE, right=FALSE)
x_SB
## class      : RasterLayer 
## dimensions : 3399, 3582, 12175218  (nrow, ncol, ncell)
## resolution : 9.316244e-05, 9.316384e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : NevTol_Pen 
## values     : 0.2, 1  (min, max)
plot(x_SB)

# save new raster
writeRaster(x_SB, filename="../spatial/Slope/rcl_S_B.tif", format="GTiff", overwrite=TRUE)
writeRaster(x_SB, filename="../spatial/Slope/rcl_S_B.asc", format="ascii", overwrite=TRUE)

#############Raster Slope Hipotesis "C"######################
rcl_S_C <- read.table("../spatial/Slope/RasterValue_SlopeC.txt",sep = ",", header=F)
dim(rcl_S_C)
## [1] 2 3
class(rcl_S_C)
## [1] "data.frame"
rcl_S_C<- as.matrix(rcl_S_C) 
class(rcl_S_C)
## [1] "matrix"
rcl_S_C
##        V1      V2  V3
## [1,] 0.00  6.5700 1.0
## [2,] 6.57 19.6999 0.2
# reclasify input raster Slope NT
x_SC<-reclassify(myrasterS,  rcl=rcl_S_C, include.lowest=FALSE, right=FALSE)
x_SC
## class      : RasterLayer 
## dimensions : 3399, 3582, 12175218  (nrow, ncol, ncell)
## resolution : 9.316244e-05, 9.316384e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : NevTol_Pen 
## values     : 0.2, 1  (min, max)
plot(x_SC)

# save new raster
writeRaster(x_SC, filename="../spatial/Slope/rcl_S_C.tif", format="GTiff", overwrite=TRUE)
writeRaster(x_SC, filename="../spatial/Slope/rcl_S_C.asc", format="ascii", overwrite=TRUE)

##############3Raster Slope Hipotesis "D"#######################
rcl_S_D <- read.table("../spatial/Slope/RasterValue_SlopeD.txt",sep = ",", header=F)
dim(rcl_S_D)
## [1] 2 3
class(rcl_S_D)
## [1] "data.frame"
rcl_S_C<- as.matrix(rcl_S_D) 
class(rcl_S_D)
## [1] "data.frame"
rcl_S_D
##       V1      V2  V3
## 1 0.0000  9.8499 1.0
## 2 9.8499 19.6999 0.2
# reclasify input raster Slope NT
x_SD<-reclassify(myrasterS,  rcl=rcl_S_D, include.lowest=FALSE, right=FALSE)
x_SD
## class      : RasterLayer 
## dimensions : 3399, 3582, 12175218  (nrow, ncol, ncell)
## resolution : 9.316244e-05, 9.316384e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : NevTol_Pen 
## values     : 0.2, 1  (min, max)
plot(x_SD)

# save new raster
writeRaster(x_SD, filename="../spatial/Slope/rcl_S_D.tif", format="GTiff", overwrite=TRUE)
writeRaster(x_SD, filename="../spatial/Slope/rcl_S_D.asc", format="ascii", overwrite=TRUE)

######################################################################


###################Vegetation Type###############################


#Raster Vegetation Type in Nevado de Toluca
myrasterVT<-raster("../spatial/VegetationType/nevado_f.tif")
myrasterVT
## class      : RasterLayer 
## dimensions : 3400, 3583, 12182200  (nrow, ncol, ncell)
## resolution : 9.313644e-05, 9.313644e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /LUSTRE/Genetica/ngalvez/Multihierarchical_NevadoToluca/spatial/VegetationType/nevado_f.tif 
## names      : nevado_f 
## values     : 0, 255  (min, max)
plot(myrasterVT)

###########Vegetation Hipotesis "A"##########
#read reclasification matrix
rclA <- read.table("../spatial/VegetationType/RasterValue_A_G.txt",sep = ",", header=F)
dim(rclA)
## [1] 7 3
class(rclA)
## [1] "data.frame"
rclA<- as.matrix(rclA) 
class(rclA)
## [1] "matrix"
rclA
##      V1 V2  V3
## [1,]  7  8 0.0
## [2,]  6  7 0.0
## [3,]  5  6 0.1
## [4,]  4  5 0.1
## [5,]  3  4 0.7
## [6,]  2  3 1.0
## [7,]  1  2 1.0
# reclasify input raster 
xVA<-reclassify(myrasterVT,  rcl=rclA, include.lowest=FALSE, right=FALSE)
xVA
## class      : RasterLayer 
## dimensions : 3400, 3583, 12182200  (nrow, ncol, ncell)
## resolution : 9.313644e-05, 9.313644e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : nevado_f 
## values     : 0, 1  (min, max)
plot(xVA)

# save new raster
writeRaster(xVA, filename="../spatial/VegetationType/rcl_A_G_3.tif", format="GTiff", overwrite=TRUE)
writeRaster(xVA, filename="../spatial/VegetationType/rcl_A_G_3.asc", format="ascii", overwrite=TRUE)


###########Vegetation Hipotesis "B"##########
#read reclasification matrix
rclB <- read.table("../spatial/VegetationType/RasterValue_B_D.txt",sep = ",", header=F)
dim(rclB)
## [1] 7 3
class(rclB)
## [1] "data.frame"
rclB<- as.matrix(rclB) 
class(rclB)
## [1] "matrix"
rclB
##      V1 V2  V3
## [1,]  7  8 0.0
## [2,]  6  7 0.0
## [3,]  5  6 0.1
## [4,]  4  5 0.1
## [5,]  3  4 0.5
## [6,]  2  3 1.0
## [7,]  1  2 1.0
# reclasify input raster 
xVB<-reclassify(myrasterVT,  rcl=rclB, include.lowest=FALSE, right=FALSE)
xVB
## class      : RasterLayer 
## dimensions : 3400, 3583, 12182200  (nrow, ncol, ncell)
## resolution : 9.313644e-05, 9.313644e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : nevado_f 
## values     : 0, 1  (min, max)
plot(xVB)

# save new raster
writeRaster(xVB, filename="../spatial/VegetationType/rcl_B_D_3.tif", format="GTiff", overwrite=TRUE)
writeRaster(xVB, filename="../spatial/VegetationType/rcl_B_D_3.asc", format="ascii", overwrite=TRUE)

###########Vegetation Hipotesis "C"##########
#read reclasification matrix

rclC <- read.table("../spatial/VegetationType/RasterValue_C_C.txt",sep = ",", header=F)
dim(rclC)
## [1] 7 3
class(rclC)
## [1] "data.frame"
rclC<- as.matrix(rclC) 
class(rclC)
## [1] "matrix"
rclC
##      V1 V2  V3
## [1,]  7  8 0.0
## [2,]  6  7 0.0
## [3,]  5  6 0.1
## [4,]  4  5 0.1
## [5,]  3  4 0.2
## [6,]  2  3 1.0
## [7,]  1  2 1.0
# reclasify input raster 
xVC<-reclassify(myrasterVT,  rcl=rclC, include.lowest=FALSE, right=FALSE)
xVC
## class      : RasterLayer 
## dimensions : 3400, 3583, 12182200  (nrow, ncol, ncell)
## resolution : 9.313644e-05, 9.313644e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : nevado_f 
## values     : 0, 1  (min, max)
plot(xVC)

# save new raster
writeRaster(xVC, filename="../spatial/VegetationType/rcl_C_C_3.tif", format="GTiff", overwrite=TRUE)
writeRaster(xVC, filename="../spatial/VegetationType/rcl_C_C_3.asc", format="ascii", overwrite=TRUE)


###########Vegetation Hipotesis "D"##########
#read reclasification matrix

rclD <- read.table("../spatial/VegetationType/RasterValue_D_A.txt",sep = ",", header=F) 
dim(rclD)
## [1] 7 3
class(rclD)
## [1] "data.frame"
rclD<- as.matrix(rclD) 
class(rclD)
## [1] "matrix"
rclD
##      V1 V2  V3
## [1,]  7  8 0.0
## [2,]  6  7 0.0
## [3,]  5  6 0.1
## [4,]  4  5 0.1
## [5,]  3  4 0.2
## [6,]  2  3 0.9
## [7,]  1  2 1.0
# reclasify input raster 
xVD<-reclassify(myrasterVT,  rcl=rclD, include.lowest=FALSE, right=FALSE)
xVD
## class      : RasterLayer 
## dimensions : 3400, 3583, 12182200  (nrow, ncol, ncell)
## resolution : 9.313644e-05, 9.313644e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : nevado_f 
## values     : 0, 1  (min, max)
plot(xVD)

# save new raster
writeRaster(xVD, filename="../spatial/VegetationType/rcl_D_A_3.tif", format="GTiff", overwrite=TRUE)
writeRaster(xVD, filename="../spatial/VegetationType/rcl_D_A_3.asc", format="ascii", overwrite=TRUE)


###########Vegetation Hipotesis "E"##########
#read reclasification matrix

rclE <- read.table("../spatial/VegetationType/RasterValue_E_B.txt",sep = ",", header=F)
dim(rclE)
## [1] 7 3
class(rclE)
## [1] "data.frame"
rclE<- as.matrix(rclE) 
class(rclE)
## [1] "matrix"
rclE
##      V1 V2  V3
## [1,]  7  8 0.0
## [2,]  6  7 0.0
## [3,]  5  6 0.1
## [4,]  4  5 0.1
## [5,]  3  4 0.2
## [6,]  2  3 0.7
## [7,]  1  2 1.0
# reclasify input raster 
xVE<-reclassify(myrasterVT,  rcl=rclE, include.lowest=FALSE, right=FALSE)
xVE
## class      : RasterLayer 
## dimensions : 3400, 3583, 12182200  (nrow, ncol, ncell)
## resolution : 9.313644e-05, 9.313644e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : nevado_f 
## values     : 0, 1  (min, max)
plot(xVE)

# save new raster
writeRaster(xVE, filename="../spatial/VegetationType/rcl_E_B_3.tif", format="GTiff", overwrite=TRUE)
writeRaster(xVE, filename="../spatial/VegetationType/rcl_E_B_3.asc", format="ascii", overwrite=TRUE)


###########Vegetation Hipotesis "F"##########
#read reclasification matrix
rclF <- read.table("../spatial/VegetationType/RasterValue_F_E.txt",sep = ",", header=F)
dim(rclF)
## [1] 7 3
class(rclF)
## [1] "data.frame"
rclF<- as.matrix(rclF) 
class(rclF)
## [1] "matrix"
rclF
##      V1 V2  V3
## [1,]  7  8 0.0
## [2,]  6  7 0.0
## [3,]  5  6 0.1
## [4,]  4  5 0.1
## [5,]  3  4 0.2
## [6,]  2  3 0.5
## [7,]  1  2 1.0
# reclasify input raster 
xVF<-reclassify(myrasterVT,  rcl=rclF, include.lowest=FALSE, right=FALSE)
xVF
## class      : RasterLayer 
## dimensions : 3400, 3583, 12182200  (nrow, ncol, ncell)
## resolution : 9.313644e-05, 9.313644e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : nevado_f 
## values     : 0, 1  (min, max)
plot(xVF)

# save new raster
writeRaster(xVF, filename="../spatial/VegetationType/rcl_F_E_3.tif", format="GTiff", overwrite=TRUE)
writeRaster(xVF, filename="../spatial/VegetationType/rcl_F_E_3.asc", format="ascii", overwrite=TRUE)

###########Vegetation Hipotesis "G"##########
#read reclasification matrix

rclG <- read.table("../spatial/VegetationType/RasterValue_G_F.txt",sep = ",", header=F)
dim(rclG)
## [1] 7 3
class(rclG)
## [1] "data.frame"
rclG<- as.matrix(rclG) 
class(rclG)
## [1] "matrix"
rclG
##      V1 V2  V3
## [1,]  7  8 0.0
## [2,]  6  7 0.0
## [3,]  5  6 0.1
## [4,]  4  5 0.1
## [5,]  3  4 0.2
## [6,]  2  3 0.2
## [7,]  1  2 1.0
# reclasify input raster 
xVG<-reclassify(myrasterVT,  rcl=rclG, include.lowest=FALSE, right=FALSE)
xVG
## class      : RasterLayer 
## dimensions : 3400, 3583, 12182200  (nrow, ncol, ncell)
## resolution : 9.313644e-05, 9.313644e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : nevado_f 
## values     : 0, 1  (min, max)
plot(xVG)

# save new raster
writeRaster(xVG, filename="../spatial/VegetationType/rcl_G_F_3.tif", format="GTiff", overwrite=TRUE)
writeRaster(xVG, filename="../spatial/VegetationType/rcl_G_F_3.asc", format="ascii", overwrite=TRUE)


###########Vegetation Hipotesis "H"##########
#read reclasification matrix

rclH <- read.table("../spatial/VegetationType/RasterValue_H_H.txt",sep = ",", header=F)
dim(rclH)
## [1] 7 3
class(rclH)
## [1] "data.frame"
rclH<- as.matrix(rclH) 
class(rclH)
## [1] "matrix"
rclH
##      V1 V2  V3
## [1,]  7  8 0.0
## [2,]  6  7 0.0
## [3,]  5  6 0.1
## [4,]  4  5 0.1
## [5,]  3  4 0.7
## [6,]  2  3 0.2
## [7,]  1  2 1.0
# reclasify input raster 
xVH<-reclassify(myrasterVT,  rcl=rclH, include.lowest=FALSE, right=FALSE)
xVH
## class      : RasterLayer 
## dimensions : 3400, 3583, 12182200  (nrow, ncol, ncell)
## resolution : 9.313644e-05, 9.313644e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : nevado_f 
## values     : 0, 1  (min, max)
plot(xVH)

# save new raster
writeRaster(xVH, filename="../spatial/VegetationType/rcl_H_H_3.tif", format="GTiff", overwrite=TRUE)
writeRaster(xVH, filename="../spatial/VegetationType/rcl_H_H_3.asc", format="ascii", overwrite=TRUE)


###########Vegetation Hipotesis "I"##########
#read reclasification matrix

rclI <- read.table("../spatial/VegetationType/RasterValue_I_I.txt",sep = ",", header=F)
dim(rclI)
## [1] 7 3
class(rclI)
## [1] "data.frame"
rclI<- as.matrix(rclI) 
class(rclI)
## [1] "matrix"
rclI
##      V1 V2  V3
## [1,]  7  8 0.0
## [2,]  6  7 0.0
## [3,]  5  6 0.1
## [4,]  4  5 0.2
## [5,]  3  4 0.2
## [6,]  2  3 0.2
## [7,]  1  2 1.0
# reclasify input raster 
xVI<-reclassify(myrasterVT,  rcl=rclI, include.lowest=FALSE, right=FALSE)
xVI
## class      : RasterLayer 
## dimensions : 3400, 3583, 12182200  (nrow, ncol, ncell)
## resolution : 9.313644e-05, 9.313644e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : nevado_f 
## values     : 0, 1  (min, max)
plot(xVI)

# save new raster
writeRaster(xVI, filename="../spatial/VegetationType/rcl_I_I_3.tif", format="GTiff", overwrite=TRUE)
writeRaster(xVI, filename="../spatial/VegetationType/rcl_I_I_3.asc", format="ascii", overwrite=TRUE)

###########Vegetation Hipotesis "J"##########
#read reclasification matrix

rclJ <- read.table("../spatial/VegetationType/RasterValue_J.txt",sep = ",", header=F)
dim(rclJ)
## [1] 7 3
class(rclJ)
## [1] "data.frame"
rclJ<- as.matrix(rclJ) 
class(rclJ)
## [1] "matrix"
rclJ
##      V1 V2  V3
## [1,]  7  8 0.0
## [2,]  6  7 0.0
## [3,]  5  6 0.1
## [4,]  4  5 0.5
## [5,]  3  4 0.7
## [6,]  2  3 1.0
## [7,]  1  2 1.0
# reclasify input raster 
xVJ<-reclassify(myrasterVT,  rcl=rclJ, include.lowest=FALSE, right=FALSE)
xVJ
## class      : RasterLayer 
## dimensions : 3400, 3583, 12182200  (nrow, ncol, ncell)
## resolution : 9.313644e-05, 9.313644e-05  (x, y)
## extent     : -99.99505, -99.66134, 18.99265, 19.30931  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : nevado_f 
## values     : 0, 1  (min, max)
plot(xVJ)

# save new raster
writeRaster(xVJ, filename="../spatial/VegetationType/rcl_I_I_3.tif", format="GTiff", overwrite=TRUE)
writeRaster(xVJ, filename="../spatial/VegetationType/rcl_I_I_3.asc", format="ascii", overwrite=TRUE)

##############################################